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We study the statistical properties of charge and energy transport in electron conducting junctions 
with electron-phonon interactions, specifically, the thermoelectric efficiency and its fluctuations. The 
system comprises donor and acceptor electronic states, representing a two-site molecule or a double 
quantum dot system. Electron transfer between metals through the two molecular sites is coupled to 
a particular vibrational mode which is taken to be either harmonic or anharmonic- a truncated (two- 
state) spectrum. Considering these models we derive the cumulant generating function in steady 
state for charge and energy transfer, correct to second-order in the electron-phonon interaction, but 
exact to all orders in the metal-molecule coupling strength. This is achieved by using the non¬ 
equilibrium Green’s function approach (harmonic mode) and a kinetic quantum master equation 
method (anharmonic mode). From the cumulant generating function we calculate the charge current 
and its noise and the large deviation function for the thermoelectric efficiency. We demonstrate that 
at large bias the charge current, differential conductance, and the current noise can identify energetic 
and structural properties of the junction. We further examine the operation of the junction as a 
thermoelectric engine and show that while the macroscopic thermoelectric efficiency is indifferent to 
the nature of the mode (harmonic or anharmonic), efficiency fluctuations do reflect this property. 


I. INTRODUCTION 

Single-molecule junctions offer a versatile playground 
for probing basic questions in condensed phases physics: 
How do quantum effects and many-body interactions 
(electron-electron, electron-phonon) control charge and 
energy transport processes, thus the operation of nano¬ 
scale atomic and molecular devices^? How do we 
accurately and efficiently simulate quantum transport 
phenomena involving different particles (and quasi¬ 
particles), electrons, phonons, spins, polarons? Recent 
progress in experimental techniques has made it possi¬ 
ble to perform sensitive measurements at the molecu¬ 
lar scale, in the linear and non-linear transport regimes, 
to observe signatures of many body effects. Kondo 
physics, the hallmark of strongly correlated electrons, 
was observed in different molecules, see e.g. Ref—. Cou¬ 
pled electron-vibration processes were probed in single¬ 
molecule junctions applying inelastic electron tunneling 
spectroscopy^— and Raman spectroscopy tools£~—, dis¬ 
playing frequency shifts and mode heating in response to 
electron conduction. Noise characteristics of the charge 
current can further expose the nature of the vibrational 
modes contributing to electron dynamics ^ 10 ' 11 . 

Theoretical and computational methodologies dedi¬ 
cated to the effects of electron-phonon interactions on 
transport in nano-conductors were reviewed in Refs: 12 ' 13 . 
The celebrated Anderson-Holstein model, with a single 
electronic orbital coupled to a local phonon mode, ex¬ 
poses an intricate interplay between the electronic and 
nuclear degrees of freedom. The model has been exten¬ 
sively studied to reveal the behavior of the current and 
its fluctuations in different regimes of electron-phonon 


coupling, see for example^ - —. An extension of the 
Anderson-Holstein model with a secondary phonon bath 
was examined in many studies, see e.g. the exploration 
of thermoelectric transport in a three-terminal junction^ 
and the analysis of transient effects— 

Complementing the Anderson-Holstein model, the 
donor-acceptor (DA) prototype junction of Fig. |T] al¬ 
lows the exploration of a broad range of problems. The 
model comprises two electronic states, referred to as 
“donor” (D) and ’’acceptor” (A), following the chem¬ 
istry literature on electron transfer reactions. Electron 
transfer between the D and A sites is assisted by a par¬ 
ticular vibrational mode, isolated, or coupled to a sec¬ 
ondary phonon bath. This model, suggested to describe 
a molecular electronic rectifier—, was recently revisited 
and analyzed using a variety of tools, for example, the 
Fermi golden rule approac h 34 : 35 , quantum master equa¬ 
tions (QME)2k2Z, the non-equilibrium Green’s function 
(NEGF) technique^— and from influence functional 
path integral simulations— In molecules, the DA model 
represents charge transfer between weakly connected 
chemical groups, facilitated by a vibrational mode— In 
the context of nanoelectromechanical systems^ the two 
electronic states can be realized by gate-controlled quan¬ 
tum dots which are coupled by a mechanical system^, a 
suspended tunnel-junction such as a carbon nanotube. 
Electrons transferred between the quantum dots can 
e.g. excite transverse acoustic modes of the suspended 
tub o 36 : 44 . The single bosonic mode can also represent a 
cavity mode assisting electron tunneling between quan¬ 
tum dots, and other hybrid model s 45 : 46 . 

The DA model in Fig. |T] offers a rich setting for inves¬ 
tigating the role of inelastic-vibrationally assisted elec- 
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FIG. 1. (color online) Scheme of a donor (D) - acceptor (A) molecular junction, (a) In the DA-AH model electron transfer is 
coupled to a highly anharmonic vibrational mode which consists of two levels, (b) In the DA-HO model the vibrational mode 
is assumed harmonic. The arrows exemplify an inelastic coupled transport process, with an electron hopping from left to right, 
assisted by an excitation of the mode. In the analysis of thermoelectric devices we set Tr > Tl, Hl > Ma¬ 


tron scattering in far-from-equilibrium (nonlinear) situa¬ 
tions. In contrast to coherent conduction, inelastic elec¬ 
tron transport can realize nontrivial effects beyond linear 
response, such as charge and thermal rectification and 
cross-rectification effects. This was e.g. demonstrated in 
Ref. 47 within a three-terminal DA junction, by replacing 
the single vibration by a phonon bath. 

Traditionally, quantities of interest in transport exper¬ 
iments and calculations were averaged values for popula¬ 
tion and currents. However, it should be recognized that 
nanoscale junctions suffer from strong random fluctua¬ 
tions due to their surrounding environments. It is there¬ 
fore desirable to develop a probabilistic theory for mea¬ 
surable quantities. Indeed, in small systems the second 
law of thermodynamics should be replaced by a universal 
symmetry, the fluctuation theorem^—, 


where Pt(S) ( Pt(—S )) is the probability distribution for 
observing a positive (negative) entropy production dur¬ 
ing the time interval t. 

Fluctuations in heat provided to a nanoscale engine 
and work performed naturally translate to stochastic ef¬ 
ficiencies. Universal characteristics of the corresponding 
probability distribution function P t (7/), for time-reversal 
symmetric engines, were explored in Refs— , from the 
principles of classical stochastic thermodynamics. It can 
be proved that the large deviation function (LDF) for effi¬ 
ciency, defined as J(rj) = — lim^oo i _1 ln[P t (77)], attains 
a global minimum which coincides with the macroscopic 
(average) efficiency, and a global maximum which corre¬ 
sponds to the least probable efficiency, coinciding with 
the Carnot efficiency^—. These universal features are 
a direct consequence of the fluctuation theorem. Exten¬ 
sions of this analysis to explore efficiency statistics for 
systems with broken time-reversal symmetry were given 
in Ref.—. Fluctuations of the finite-time efficiency were 
experimentally demonstrated in Ref—. Beyond classical 
thermodynamics, in a recent study the concept of the 
“stochastic efficiency” was examined within a quantum 


coherent model of a thermoelectric junction, by employ¬ 
ing the non-equilibrium Green’s function technique- 

in this paper, we provide a complete analysis of charge 
and energy transport behavior in the donor-acceptor 
junction by taking a full counting statistics (FCS) ap¬ 
proach. We consider two variants of the model as de¬ 
picted in Fig. [I] Electron transfer may couple to a har¬ 
monic vibrational mode (DA-HO model), or to an anhar¬ 
monic impurity (DA-AH model). The latter case is rep¬ 
resented by a two-state system, a truncated vibrational 
manifold. 

We rigorously derive the cumulant generating func¬ 
tions for the models of Fig. [I] under the assumption of 
weak electron-vibration coupling, handing over the com¬ 
plete information over the models’ steady state trans¬ 
port behavior. From the cumulant generating functions 
we explore transport in the junctions far from equilib¬ 
rium, specifically, we aim in identifying transport quan¬ 
tities which are sensitive to the nature of the vibrational 
mode. Finally, we investigate thermoelectric efficiency 
fluctuations in the DA junction beyond linear response. 
Interestingly, our analysis in this paper exemplifies that 
one can reconcile two central yet disparate techniques, 
QME and NEGF, to obtain consistent results, within the 
same order in perturbation theory. 

The paper in organized as follows. We introduce the 
DA junction in Section |TT] and perform a FCS analysis 
in Sec. mu The cumulant generating function (CGF) of 
the DA-AH model is derived in Sec. IIII Al by employ¬ 
ing the quantum master equation approach. The deriva¬ 
tion of the CGF for the DA-HO case is detailed in Sec. 
IIII B1 using the non-equilibrium Green’s function tech¬ 
nique. Two applications are described in Sec. HVl We 
simulate the junction’s current-voltage characteristics in 
Sec. IIV Al and examine the statistics of the thermoelec¬ 
tric efficiency in Sec. IIV Bl further comparing numeri¬ 
cal results far from equilibrium with the linear response 
(Gaussian) limit. Our findings are summarized in Sec. 
m For simplicity, we set e = h = kg = 1 throughout 
derivations. 
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II. MODEL 

The DA junction includes a two-site structure, rep¬ 
resenting a donor-acceptor molecule (or equivalently, 
a double-quantum-dot system), placed in between two 
metal leads. The total Hamiltonian is given by 

H t = H m + H l + H r + H c + H vib + H f . (2) 

The molecular Hamiltonian Hm includes the donor and 
acceptor sites, 

H m = e d c\c d + e a ctc a , ( 3 ) 

with Cd/ a (c^ a ) as a fermionic annihilation (creation) op¬ 
erator at the donor or acceptor sites with energies e d / a - 
The second and third terms in Eq. © represent the left 
(Hl) and right (Hr) metal leads, modeled by collections 
of non-interacting electrons, 

H l e/CjCj, -Hr = ^ c r c\.c r , (4) 

leL r£R 

with the fermionic annihilation (creation) operators Cj 
(c'j). The tunneling energies of electrons from the donor 
(acceptor) site to the left (right) lead vi (v r ) are included 
in He , and are assumed to be real valued, 

He = y; vijcjcd + c\ci) + y t v r (cic a +4c r ). (5) 

leL r£R 

H v ib and Hj represent the Hamiltonians of the molecular 
vibrational mode and its coupling with the D and A sites 
(strength g). Assuming an harmonic local mode we write 

H v ib — 

Hi = g[c\c a +clcd\(bl + bo), (6) 

with & 0 (&J) a bosonic annihilation (creation) operator for 
a vibrational mode of frequency u> o. The interaction Hi is 
sometimes refereed to as an “off-diagonal” since electron 
exchange between the two sites is allowed only via the ex¬ 
citation and/or relaxation of the mode, see Fig. [Q Note 
that we do not include here a direct tunneling (elastic) 
term between sites D and A. This simplification allows 
us to reach closed analytic results for the CGF. The con¬ 
tributing of elastic tunneling processes can be included 
in an additive manner—, a reasonable approximation at 
weak coupling— 

We diagonalize the electronic Hamiltonian, H e \ = 
Hm + Hl + Hr + He, and write it down in terms of 
a new set of fermionic operators (ai/ r , aj^), 

H e i = y^iajai + y ^e r al.a r . (7) 

l r 

These operators are related to the original set by2Z 
Q = y^7 m, c a = r a r , 

l r 

Cl = ^2 Cr = ^2 ^rr'dr', (8) 

l' r' 


with the dimensionless coefficients 


7 1 = 


vi 


£; Cd e, — 


mi' = 8u> - 


vav 


V ei— e t i+iS 


ei - e v + 


iS’ 


(9) 

Here S is a positive infinitesimal number introduced to 
ensure causality. Analogous expressions hold for the r 
set. The expectation values of number operators obey, 
e.g., at the L end, (a\av) = SwfL(ei) with / l(ci) = 


{ exp[/3i(e; — hl)] + l} 1 the Fermi distribution for the 
lead with chemical potential /jll and an inverse temper¬ 
ature Pl- Using the new operators, we write down the 
total Hamiltonian for the DA-HO junction as 


Hda-ho = H e i + u} 0 blb 0 

+ 9 \lllr^[a r + h.c.] (&o + &o)-(10) 

iei.rei? 


The last term in the Hamiltonian describes electron-hole 
pair generation assisted by the interaction with the vi¬ 
brational mode. 

In a simpler version of this model we replace the infinite 
level spectrum of the harmonic oscillator by a truncated 
two-level system, to mimic a highly anharmonic vibra¬ 
tional mode. The Hamiltonian of this DA-AH model can 
be conveniently written in terms of the Pauli matrices as 


H D a-ah = H e i + —~(Tz + g ^2 [lilra\a r + h.c.]a x . 

l£L,r£R 

( 11 ) 

The DA-HO and DA-AH models described above are 
simple enough to allow us to derive expressions for the 
corresponding CGFs. Meanwhile, we (i) reconcile the 
NEGF approach with QME, to assist in method devel¬ 
opment in the area of quantum transport, and (ii) provide 
analytic results for transport characteristics, to bring in¬ 
tuitive guidelines for functionality. 

We substantiate our modeling by describing connec¬ 
tions to ab-initio studies of molecular conduction— The 
donor and acceptor sites could represent spatially sep¬ 
arated electronic states in a molecule or atoms in an 
atomic wire, with the electron-phonon interaction of Eq. 
© describing a bond-length stretching mode. The rigid 
motion of a molecule/atomic chain between the metal 
leads can be modeled by an additional electron-phonon 
interaction Hamiltonian H^ = gi(b\ + 6i)(c^c<i + c^ a c a ). 
This mode is not included in the present study. As well, 
we ignore direct tunneling between electronic sites in 
the form H tun n = Vd, a (c^Cg + clcd^j • Detailed ab-initio 

studies of electron-phonon inelastic effects in molecular 
junctions prepare direct tunneling elements, frequencies 
of active modes, and their electron-phonon matrix ele¬ 
ments, see e.g. Refs— r— . Certain type of molecular 
junctions could be well represented by our model, when 
electron transport due to Hi dominates over both elastic 
tunneling H tU nn and the diagonal electron-phonon cou¬ 
pling Hj 1) . This is the case e.g. in Ref—, consider¬ 
ing charge transfer through a biphenyl molecule with the 
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torsion motion assisting electron hopping between the 
(almost orthogonal) benzene rings. Particularly, when 
ei / e a , calculations of transport in double quantum dot 
systems show that the inelastic component of the current 
can dominate the elastic term^i. 

We also justify our model in the language of molecular 
orbitals, electronic eigenstates of the molecule. Consider 
two orbitals, each coupled to both metal leads- but in 
an asymmetric manner: one orbital couples strongly to 
the left lead but weakly to the right, the other molecular 
orbital is strongly coupled to the right lead but weakly 
to the left side. In the absence of electron-phonon inter¬ 
action this molecule supports very small currents. It will 
however turn into a good conductor at high enough tem¬ 
peratures when phonons contributing to (@ are active, 
supporting inelastic current. For an extended discussion, 
see Ref^. 

Turning to the the AH model, the two-state impurity 
describe deviations from the harmonic picture. Besides 
electron-phonon coupled situations, the model could rep¬ 
resent electron transport junctions interacting with a lo¬ 
cal spin impurity, see e.g. Refs . 64 ! 65 , demonstrating elec¬ 
tronic read-out of nuclear spins. 

We finally comment that in the non-crossing approx¬ 
imation, when elastic and inelastic tunnelling events do 
not interfer e) 39 ! 66 , the current (or more generally, the cu- 
mulant generating function) can be written as a sum of 
elastic and inelastic terms. While here we treat the in¬ 
elastic component only, the CGF for elastic transport is 
well known§2, bringing in the standard Landauer formula 
for currents. 

In the next Section we develop a counting statistics 
approach for charge and energy transfer processes. We 
then analyze first the (simpler) DA-AH model using a 
QME approach, followed by the investigation of the DA- 
HO junction by utilizing the NEGF technique. 


III. COUNTING STATISTICS FOR CHARGE 
AND ENERGY CURRENTS 


The cumulant generating function contains informa¬ 
tion over the statistics of transferred particles and energy 
flowing across the system, potentially far from equilib¬ 
rium. Here we are interested in the CGF for coupled par¬ 
ticle (charge) and energy currents. Such a two-parameter 
CGF is necessary for obtaining later in Sec. IIV Bl the 
statistics of efficiency in a thermoelectric engine. 

We define the particle (p) and energy (e) current op¬ 
erators from the rate of change of electron number and 
electron energy in one of the leads, say R , and write 


I P (t) = ~ 


dNg(t) 


dt 


Ie(t) = - 


dt 


( 12 ) 


Here N R = E reR 4«r is the number operator for the 
total charge in the right compartment (right lead plus 
attached acceptor site). Similarly, H R = Erefl er '4 a »~' 


The operators are written in the Heisenberg (H) pic¬ 
ture, thus they should be evolved with the total Hamil¬ 
tonian for either model, A H (t) = U^(t) AU(t) where 
U(t) = e~ lHTt . We follow the convention that the cur¬ 
rent flowing out of the right lead is positive. Changes 
in the total energy and electron number in the R lead 
during the time interval (to = 0, t) (to and t are initial 
and final observation time, respectively), are given by the 
integrated currents. 


Qe(Mo) = f h(t') 
Jt o —0 

Q P (t,t 0 )= f i P (t') 

J £n —0 


dt' 


dt' 


H R (0)-H*(t), 
N R (0)-Ni(t). (13) 


Since at any instant [N r ,H r ] = 0, it is possi¬ 
ble to construct the so-called “characteristic function” 
Z(X e ,X p ), corresponding to the joint probability distri¬ 
bution P(Q e , Q p ) for the charge and energy currents. 
Following the two-time measurement procedure one can 
define the characteristic function a s 49 i 50 i 68 i 69 


Z(X e , X p ) = ^ e iXeHR+iXpNR e~ iXeH * ( t '>~ iX P N R W^, (14) 


where X e and X p are the counting fields for energy and 
particles, respectively. (• • •) represents an average with 
respect to the total density matrix pt(0) at the initial 
time. We assume that p R ( 0 ) = Pl( 0 ) < 8 > Pr(0) < 8 > Pvib(0), 
a factorized-product form for the electronic degrees of 
freedom and for the vibrational part. The leads are 
maintained in equilibrium at temperature T a = 1 /fd a 
and chemical potential it a , a = L,R, and the states 
are described by the grand canonical distribution func¬ 
tion p a (0) = exp[-/3 a (H a - p a N a )]/Z a , with Z a = 
Tr[exp[— fd a (H a — p a Na)}] as the grand canonical parti¬ 
tion function. Eq. m can be reorganized as 

Z(A e , X p ) = Tr T [t/_ Ae/2 ,_ Ap/2 (t) pr( 0) ^. /2 , Aj(/2 (*)], 

= Tl 'T [Px et x p {tj\, 

= Tr vib [pt,xM- ( 15 ) 

The second line introduces the definition of the total, 
counting-field dependent, density operator. We trace out 
its electronic degrees of freedom, (Tr e i) and express the 
characteristic function in terms of the reduced density 
matrix p A * h Ap ( t ) for the vibrational mode, 

pf» Xp (t) = Tr el [[/_ Ae/2 ,_ Ap/2 (f) M0) Ul /2iXp/2 (t )\. 

(16) 


Note that the forward and backward evolution operators 
are not hermitian conjugates. For example, the forward 
propagator is 


U- Ae/2,—Ap/2W — 


exp 


Nr 

2 2 


U(t) exp 


i^Hr + ^Nr 


= exp HH_ Ae/2 ,_ Ap/2 (f)], 


(17) 
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with the counting-field dependent total Hamiltonian 


II 


— A e /2, —Ap/2 


exp 


-i^-H R - i^-N r 
2 2 


Ht exp 


i^Hn + i^Nn 


= H e \ + H vib + S <g> [ 5 ^ 7 ; * 7 r a / t a r e 2 ( Ae£r+A p) + h.c.]. 

l,r 

(18) 


S' is a system operator; the Hr a-ah model is reached 
when S = a x , H v n, = 7 ^ 07 . The model Hr a-ho is real¬ 
ized with S = (&J + & 0 ) and H vR) = wohj&o- The electron- 
phonon coupling term here depends on the counting field. 
Herein, we use A as a short-hand notation for both A e 
and A p . To facilitate our discussion below, we define the 
operators B± x / 2 as 


B^x /2 = g[J2'yi'y r a]a r e ± ^ £rK+Xr ‘ ) +h.c.]. (19) 

l,r 


These operators correspond to the bath operator coupled 
to the system, see Eqs. m and ED- now dressed by 
counting fields A e ,A p as a consequence of the measure¬ 
ments of charge and energy. Note that the sign conven¬ 
tion for B corresponds to the respective time evolution 
operator. The counting-field Hamiltonians [Eq. (fl8l) and 
the complementarity term for the backward evolution] 
can be organized in a form convenient for a perturbation 
expansion in g , 


H± A/2 — H 0 + V±A/2> 

Hq = H e \ + H V i b , V±a /2 = S ( 8 > B±\/ 2l (20) 

I 


with H e 1 incorporating the two metals with the hy¬ 
bridized states, see Eq. 0. 


A. DA-AH model: Quantum master equation 
approach 


We derive the counting field dependent quantum mas¬ 
ter equatio n 49 ' 71 for the model ED under the assumption 
that the coupling between electron-hole pair generation 
and the vibrational mode is weah£2. Unlike other stan¬ 
dard QME approaches for molecular junctions, in which 
the molecule-metal coupling is considered wea k 56 ! 73 , in 
the present derivation the metal-molecule hybridization 
[defined below Eq. (l27l) ] can be made arbitrarily large, 
absorbed into the leads spectral density by the exact di- 
agonalization procedure presented in Sec. [Til Taking the 
time-derivative of Eq. ED we get 


P\ b (t ) = Trei [-iV-\/2{t)pl(t) + ipl(t)V x/2 (t)\ .(21) 


The operators here are written in the interaction repre¬ 
sentation, A(t) = e lHot Ae ~ lHot , with Hq including the 
uncoupled electrons and vibration. By formally integrat¬ 
ing this equation we receive the exact form {{B± x / 2 ) = 
0 ), 


P v x(t) = ~ 


dt'Tie 


to—0 


V-\/ 2 (t)V-\/ 2 (t')P\ (O + Px {t')V\/ 2 (t')V\/ 2 {t) 


-V-\/ 2 (t')p x (t')V x/2 {t) - V_ x/2 {t)p x (t')V x/2 {t') 


( 22 ) 


We now follow standard steps as in the derivation of 
the weak coupling-Markov Redfield equation. The initial 
condition is assumed fully factorized, p x ( t' ) is replaced by 
the initial condition, pj(0) = pt(0), and the upper limit 
of integration in extended to infinity, assuming Marko- 
vianity of the electron baths. The equation of motion for 
the counting field dependent reduced density matrix (de¬ 
scribing the dynamics of the vibrational mode) depends 


on the following relaxation (kd) and excitation ( k u ) rates 

/ OO 

dre-“ oT ( J BA/2(0)H_A/ 2 (T)) el , 

-OO 

/ (X) 

dre iuioT (B x/2 (0)B_ x/2 (T)) el 

-OO 

= kd[oj 0 —> —wo]- (23) 

Here (• • • } e i = Tr^Tr^- ■ • pl(0)pr( 0)]. An explicit cal¬ 
culation of the relaxation rate gives 


kd = 2ng 2 [^2\7i\ 2 \7r\ 2 fL{ei)(l ~ /fl(er))e l(Ap+erAe) ^(ei - e r +w 0 ) 

l,r 

+ M 2 |7r| 2 /fl(er)(l - /i(e;))e l(Ap+£rAe) 'T(ei ~ e r - w 0 ) . 

l,r 


( 24 ) 
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The first term represents an inelastic process with an 
electron hoping from the left lead to the right one, by 
absorbing one quanta too, satisfying energy conservation 
with e r = ei + ujq. This process, which goes against 
our convention of a positive current (flowing right to 
left), contributes negative charge and energy currents, 
reflected by the negative sign for A p and X e in the ex¬ 
ponent. The second term in Eq. (l24l) corresponds to the 
reversed process with electron hopping from the right 
metal to the left, observing e; = e r + ui o, with a positive 
sign for X p and X e . This downward rate assists in cooling 
the vibrational mode. In the complementary excitation 
process electrons lose energy to the vibration, heating 
up the junction. Eq. (l24l) can be decomposed into two 
separate contributions, 

k X d = [k X d ] L ^ R +[k X d } R ^ L , (25) 

and similarly for fc A . We define the spectral densities for 
the baths (metals) as 


J a (e) = 2 tt 5 ^ | 7 fc | 2 5 (e - e fe ). (26) 


Using the transformations © we note that these func¬ 
tions are Lorentzian-shaped, centered around the donor 
(ed) or acceptor (e a ) site energies, 


7 7 ^ ^(e) 

L[) 5 (e-e d ) 2 +r L (e) 2 /4 

T / \ _ r fl( £ ) _ 

* () 9 (e-e a )*+T R (e) 2 /4 


(27) 


with T a (e) = 27r^) fc6a n 2 <5(e — < 7 ). In terms of these 
spectral functions, bath induced relaxation rates (1551) are 
given by 


[$] L ^ R = j ^[fL(e)(l-f R (e + uJo))Me)Me + u;o) 


—i(A p -(-Ce-}-£j;o) A e 


xe 
1 A] 


[K 


xe 


de 

2i t L' 


-(Ap+eAe) 


/fl(e)(l ^ /i(e +uj 0 ))J R (e)J L (e + u> 0 ) 

(28) 


We also define the A = 0 rates from Eqs. (1551) - (1551) . only 
missing the A identifier. The rates are nonzero as long 
as: (i) for a given energy, both left and right leads are 
not fully occupied or empty and (ii) the overlap between 
the spectral densities differing by one quanta of energy is 
non-negligible. Note that because of the weak electron- 
phonon coupling approximation, each electron tunnelling 
process involves absorption/mission of a single quanta of 
energy. In other words, the dynamics is completely de¬ 
scribed by single-phonon excitation and relaxation rates. 

As mentioned above, we apply the weak-coupling Born 
Markov approximation on Eq. m- By further ignoring 
off-diagonal coherence elements for the reduced density 


matrix, we obtain the population dynamics for the vi¬ 
brational states (written here for an arbitrary number of 
levels, m = 0,1, 2 ,...), 

Pm(t) = - [mk d + (m + 1 )k u ] p^(t) 

+ (to + l)fcd pA+iM + mk uPt- iW> (29) 

where p^(t) = (m\p v x lb (t)\m) and |to) denotes the TO-th 
vibrational level. k d and k u are rates evaluated at A = 0. 
For the DA-AH model, to = 0,1, it can be shown that off- 
diagonal coherences do not appear in the Born-Markov 
approximation without further assumptions!!, and the 
dynamics of the population follows 

PoO) = - k uPo(t) +^Pi (4 

Pi(t) = kuPo(t)-k d Pi(t). (30) 

These equations can be written in a matrix form as 


^M = £(A)|/(t)), (31) 

where |p A (f)) = (Po(t),PiW)- The long-time (steady 
state) limit defined as 

G{ A) = lim ^lnZ(A) = lim - ln(/|p A (t)), (32) 

provides the CGF, where (I| = (1,1) T is the identity 
vector. In this limit, only the smallest eigenvalue of the 
Liouvillian survives. The final result for the CGF for the 
DA-AH junction is given by 

Qah (A) = — — (feu + k d ) + — \J (k u — k d ) 2 +4 fc A fc^. 

(33) 

We label this CGF by ‘AH’, to highlight the mode an- 
harmonicity. Recall that A collects two counting fields 
\ p and A e , for charge and energy, respectively. The CGF 
satisfies the fluctuation symmetry 


Gah (A e , A p ) 

= Gah(— A e + z(/3z, — Pr), —Xp + i(/3 R p R — Plpl)), 

(34) 


which can be verified by examining the rates in Eq. (1241) . 
Under the transformations A e — >• —A e + i{Ph — Pr) and 
X p —> — X p + i{P R p R — PlI-Ll) the rates transform as 

k X d -»■ fc A e^"°, 

fc A ^fc A e-^“°. (35) 

The extra factors e ±lJLUl0 cancel out in Eq. (1531) . to satisfy 
the fluctuation symmetry. 

The charge and energy currents and the corresponding 
zero frequency noise powers can be readily obtained by 
taking derivatives of the CGF with respect to the count¬ 
ing fields. For example, the particle (p) and energy (e) 
currents are obtained from 


{I P ) = 

(Ie) = 


( Qp) 

t 

(Qe) 

t 


dG(X e ,X p ) 

d(iX p ) 


A e —Ap —0 


dg(A e ,Ap) 

d(iX e ) 


A e —Ap —0 


(36) 
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The zero frequency noise of these currents are 


(S P ) = 


(( Q 2 P )) 

t 


{Se) = 


((QD) 

t 


d 2 G(X ei X p ) 
d(i\ P ) 2 
d 2 G(X e ,X p ) 
d(iX e ) 2 



(37) 


where (( Q 2 tP )) = {Qt, P )~(Qe, P ) 2 is the second-cumulant. 


B. DA-HO model: Non-equilibrium Green’s 
function approach 


We now focus our attention to the DA-HO junction of 
Fig. mb), described by Eq. (flOl) . In this case, electron- 
hole pair excitation is coupled to a harmonic oscillator 
mode. We employ the NEGF techniqu e 74 ’ 75 to derive the 
CGF of this model, again correct up to the second order 
of the electron-phonon coupling parameter g. Note that 
the presence of an infinite number of vibrational levels 
for the HO mode makes it difficult to obtain a closed 
form for the CGF under the QME approach, since the 
Liouvillian is an infinite-dimensional matrix. 

We begin with Eq. (1151) , now identifying the initial time 
by foj and write down the characteristic function on the 
Keldysh contour (see Fig. [2]) as^ 


Z(X e ,X p ) — (^ Xe / 2tXv / 2 {t)U_ Xe / 2 _ Xp / 2 {t)^ 
= Tr|p T (0) T c e _i ^. 


The forward upper (backward lower) branch of the 
Keldysh contour corresponds to the modified unitary 
evolution U_ Xe/2 _ Xp/2 (t) {Ul J2 Xp/2 (t)). Evolving the 
counting fields on two branches of the Keldysh con¬ 
tour with two different signs is the main essence of 
counting statistics problems. The normalization con¬ 
dition is trivially satisfied with Z( 0,0) = 1. In the 
above expression T c is the contour ordered operator, 
which orders operators according to their contour time 
argument; earliest-time operators appear at the right. 
A (t) = (A e (r), A p (r)) is the contour time dependent 
counting function. In the upper (+) branch, A + (t) = 
(A+(t), A+(f)) = (—A e /2, —Ap/2), in the lower (-) 
branch, A ~(t) = (A“(f), A“(f)) = (A e /2,A p /2). Moving 
to the interaction picture with respect to the Hamiltonian 
Hq = H e i + H v ib , we write the characteristic function as 


Z(A e ,A P H(T c exp(-i [ drK A(T) (r))), (39) 

J C 

where the counting field dependent interaction term is 
I/ A(T) (r) = [Mt) + &J(t)] 

x g Y. 7*7r e~ t( - Xp( - T ^ +erXe ^ a] (r) a r (T) + h.c. 

l,r 

= 9^2,1*7r aj(T)a r (T) [b 0 (r) + fej(r)] + h.c. 

l, r 


Tl -A e /2, -Xp/2 

- - 


T2 A e /2, Ap/2 

FIG. 2. Keldysh contour representing the counting statis¬ 
tics problem. The forward upper (+) and backward lower (-) 
branches evolve with different Hamiltonians corresponding to 
different counting fields, n, t 2 are the contour times and t is 
the final observation time. 



introducing the short notation d r (r) = 

e~»(A p (T)+e r A e ( T ))a r (t). Our objective is to calculate the 
CGF, correct up to second-order in electron-phonon 
coupling but exact to all orders in the metal-molecule 
hybridization. It can be shown that a naive perturbative 
calculation of Eq. (l39l) in terms of the electron-phonon 
coupling g leads to a violation of the non-equilibrium 
fluctuation symmetry. Moreover, such a perturbative 
treatment cannot capture the correct non-equilibrium 
phonon distribution, and it will lead to (an incorrect) 
long-time solution for the vibrational density matrix 
which depends on the initial- arbitrary state p V ib( 0). 

In order to restore the fluctuation symmetry, and 
obtain the correct non-equilibrium phonon distribution 
[which depends on the temperatures of the electronic 
baths and their chemical potentials, but not on p n ,ib(0)], 
one has to sum over an infinite subclass of diagrams in 
this perturbative expansion. This procedure takes into 
account all electron scattering processes which are facili¬ 
tated by the absorption or emission of a single quanta wo ■ 
Physically, the summation collects not only sequentially 
tunneling electrons, but all coordinated multi-tunneling 
processes, albeit with each electron interacting with the 
mode to the lowest order, to absorb/emit each a single 
quanta u)q. This summation can be achieved by exploit¬ 
ing the random-phase approximation (RPA) as done in 
RefSi 27 ’ 76 ’ 77 , also referred to as the self-consistent Born 
approximation^. Summing over a particular set of dia¬ 
grams (ring type) in the perturbative series, see Fig. [3] 
we reach the following expression, 

lnZ RPA (A e , A p ) = — ^ Tr T In [i - D 0 (t, t')F(t', r)] . 

(40) 

The symbol Tr T denotes an integration over contour time 
variables (r, r'). For example, 

Tr T [Do(t, ,t )] = J dr J dr'D^r, t')F(t', t). 

(41) 

Here, / is the identity matrix in the Keldysh space and 
Dq(t,t') is the free phonon Green’s function, 

A>(ti,t 2 ) = —*(T c X(ti)X(t 2 )), 


(42) 














with X = (bo + 6g), proportional to the phonon displace¬ 
ment operator. F(t,t') is the counting-field dependent 
electron-hole Green’s function. It describes electron hop¬ 
ping processes from the left to the right lead, and vice 
versa, 

F(n,T 2 ) =-ig 2 M 2 7r-| 2 

leL,reR 

X [9i(Ti,T 2 )g r (T 2 ,T 1 ) + g r (T 1 ,T 2 )gi(T 2 ,T 1 )]. 

(43) 

This expression is symmetric under the exchange of the 


contour time parameters t\ and t 2 . Recall that the tilde 
symbol advices that the Green’s function is A dependent. 
This propagator involves free electron Green’s functions 
for the left and right leads, 

9i{ti,t 2 ) = —i {T c ai(r{)a\ (t 2 )), 

9r(ri,T 2 ) = -i{T c a r (n)al(T 2 )) . (44) 

Explicit expressions for different components of these 
Green’s functions in real time are given in Appendix A. 
Here, we write down the lesser component, given as, 


F K (u) = —i 2 tt^ 2 l7z| 2 |7r| 2 /i(ei)( 1 ~ /fl(er))e - e r - w) 

l,r 

+ E h\ 2 hr\ 2 f R (e r )(l - f L (e l ))e i ^ + ^S(e l e r + uj . 

l, r 


(45) 


The greater component is obtained from F > (uj) = 
F < (—uj). It is clear that the greater (lesser) compo¬ 
nent corresponds to the electronic bath-induced transi¬ 
tion rates within the vibrational mode, ( k x ), see def¬ 
initions (E3l) and a physical explanation below Eq. (1241) . 


with a counting field independent term D 0 1 . The matrix 
D A ’ 1 can be written explicitly as 


£>a>) = 


_([di ;]-»-*» 


F> M 


F<( W ) 

-[D^i^-F^u) 



(■ b ) 


Q 

y 



where stand for the retarded, advanced, 

time-ordered, anti-time ordered, lesser and greater com¬ 
ponents of the Green’s functions, respectively. The free- 
phonon Green’s functions D o’ a (w) are given by 

floM= r y?2 - 2. floM = PoM]*, (47) 

(w + ir)Y - lo q 



FIG. 3. Ring type Feynman diagrams in contour time, (a) 
Second-order, (b) fourth-order, and (c) sixth-order diagrams 
in the electron-phonon coupling. The dotted line repre¬ 
sents the phonon Green’s function Do. Closed loops are the 
electron-hole propagator F(ti,t 2 ), the sum of two diagrams 
(d) consisting of the bare left (solid) and right (dashed) leads 
Green’s functions. 

Projecting Eq. (141)1) to the real time and invoking the 
steady state limit by taking to —> — oo, we write the CGF 

a a 68.75 

£(A e ,A p ) = lim -lnZ R p A (Ae,Ap) 

t—±oo t 

= ~ j ^lndet (46) 

where D x (uj) = D^ 1 ^) — a z F(uj)cr z with a, as the 
third Pauli matrix. Note that we renormalized the CGF 


where r\ is an infinitesimal positive number, introduced 
to preserve the causality of the retarded Green’s func¬ 
tion. The determinant of D 7 (a;) can be immediately 
evaluated, 


det D x 1 (w) 



A a (w). 

(48) 


Here 


A x (u) = F<{u)F>(u) 


(f<h+f>h ) 2 

4 


(49) 


is written solely in terms of the electron-hole Green’s 
function, and it will end up being the central quantity 
in this problem. The above determinant was simplified 
using the identity F*(w) + F* (to) = F > (ou) + F < (oj ), see 
Appendix A. We now take the A derivative of the CGF 
in Eq. (l46l) . to obtain 


dxG(X) = - 


du 


d\A\(uj) 


47T 


2ujq 


(AA)] t4H 


(50) 
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The integration can be performed to the leading order in 
the electron-phonon coupling g. To the lowest nontrivial 
order in the electron-phonon coupling the location of the 
poles can be approximated by 

± |w 0 + Re[F r (a;o)] ± iyj Ax(u 0 )} , (51) 

where we identify Re[E r '] = ( F 4 — F*)/2. Employing 
the residue theorem, the integration in Eq. ED can be 
performed, resulting in 

dxG(X) « -dxx/M^j, (52) 

We now formally identify the lesser and greater compo¬ 
nents of the electron-hole Green’s functions in Aa(wo) 
with the excitation and relaxation rates, defined in the 
QME approach in Sec. IIII A1 see Eq. (HM1) . They are given 
by F > ( wo) = —ikd and F < (u)o) = —ik u . In the presence 
of counting fields, F > (uj 0 ) — — ikF < (u> 0 ) = —ik *. The 
final expression for the CGF is 

Gho(X) = ^{k d - k u ) - ^(k u + k d ) 2 - 4 k*k%. 

(53) 

Remarkably, the expressions for GHo(X e , X p ) and 
Gah(X b , X p ) in Eq. iTTTTTl) are very similar, besides the 
sign differences, though they were obtained via two com¬ 
pletely different approaches. The sign difference reflects 
the different normalization, with the AH mode conserv¬ 
ing population in the two states of the mode, while in the 
harmonic mode all levels are occupied at nonzero temper¬ 
ature. The RPA approximation restores the fluctuation 
symmetry, Eq. m- 

The CGFs for the DA-AH and DA-HO models, Eqs. 
(l3d]l and (l5dl) . respectively, are the main analytical re¬ 
sults of this paper—. It should be emphasized that only 
few impurity models, essentially, variants of the single 
dot Anderson modeU^— — , can be solved analytically, 
within certain approximations, to provide the CGF and 
expose charge and energy statistics under interactions. 
Our work here substantially extends these efforts, by 
solving the FCS of a vibrationally assisted two-site elec¬ 
tronic conduction. In Sec. IIV Al we derive and simulate 
charge and energy currents and their noise using Eqs. 
(l36l) and (IdTl) . As a further nontrivial application, we 
employ the CGFs in Sec. HVBI to simulate fluctuations 
of the thermoelectric efficiency. 


IV. APPLICATIONS 

A. Charge current and Fano factor far from 
equilibrium 

We present here simulations for the charge current 
and its noise in the DA junction. Particularly, we ex¬ 
amine signatures of mode harmonicity in transport. In 
Ref— we further explore fingerprints of vibrational an- 
harmonicities in linear response quantities: the electronic 


thermal conductance, the thermopower, and the thermo¬ 
electric figure of merit. 

We obtain closed-form expressions for the currents and 
high order cumulants by taking partial derivatives of the 
CGF with respect to the counting fields, see Eqs. m 
and ED- The particle (p) and energy (e) currents in the 
DA-AH and DA-HO models are given by 

i uR — yLuR — yL h,L — yRuL — yR 

(^jAH/HO^ _ 9 h d h u _ h d h u 

p k d + sk u 

^jAH/HO'j _ [d(i\e)ku U=o] + k u [*%A e ) k d | A=o] 

k d + sk u 

(54) 


where s = +1 (s = —1) for the AH (HO) mode. Note 
that we did not simplify the expression for the energy cur¬ 
rent above; the derivatives return energy transfer rates 
which are analogues to Eq. (1251) . only with an addi¬ 
tional energy variable in the integrand. The average heat 
current, extracted from the right terminal, is defined as 
{I^ H ^ HO ) = (lf H ^ HO ) — idR(Ip H ^ H °). We further write 
down closed expressions for the particle-current noise, 


{s ah/ho } 


2 s 

k d + sk u 


{ IA H /H°)2 


k d + ski 


( b.L—YRb.L—YR fR—yLfR—yL\ 
\^u ^d ' ^d ) • 


(55) 


The first term hands over a strictly non-equilibrium 
(finite-bias) noise. The second term survives even when 
the bias voltage is zero, thus we refer to it as the equilib¬ 
rium contribution (though it is somewhat modified with 
bias). At low bias and low temperatures, k u -C k d , thus 
{Sp H ) ~ (Sp°). At finite bias, significant differences 
show up, as we discuss in the text following Figs. [ 6 ] and 
□ 

The Fano factor, defined as the ratio of the noise to 
the current, F = ( S p )/(I p ), receives a rather transparent 
form 


pAH/HO 


+ 2 


2s(Ip H ' HO ) 
k d + sk u 

J~R — yLl,R — yL i h,L — yRuL 

rvj ft,,. ~r /tj a,,. 


R 


l,R—yLl,R—yL _ uL—yRuL—yR 

K'd ^u K'd Ni 


(56) 


The second term here does not depend on the mode har- 
monicity/anharmonicity: At finite bias k R ^ L > k^ R , 
thus this term roughly takes on the value 2 , besides at 
asymptotically small biases when the denominator drops 
to zero since the current itself is diminishing. The first 
term in Eq. (1561) . in contrast, depends on the mode har¬ 
monicity, it strongly varies with the bias voltage, and 
physically it corresponds to the ratio between two rates: 
charge transfer through the junction and transitions be¬ 
tween vibrational states within the attached mode. 

Before presenting results at finite voltage and temper¬ 
ature, we derive scaling relations for the current and 
its noise in the large voltage - zero temperature limit, 
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A g [eV] Ag [eV] 




A g [eV] Ag [eV] 


FIG. 4. (color online) Charge current and differential conductance as a function of voltage bias for the DA-AH model with 
(a)-(b) cjo = 0.1 eV, and (c)-(d) u>o = 0.05 eV. Other parameters are e d = ta = 0.25 eV, T = Tl = Tr = 100 K, g = 0.1 eV, 
T =0.01 eV (solid), T = 0.1 eV (dashed). 


when metal-molecule hybridization is large, Jl(w) —> j=- ; 
Jr(oj) —> see Eq. (E71) . This limit will allow us 

to pinpoint on fundamental differences between the HO 
and the AH mode models. As well, scaling relations will 
be contrasted with results from the Anderson-Holstein 
model. We introduce the notation g 2 = 8g 2 /tt, assume 
zero electronic temperature and a large voltage bias, 
Ag = g R - g L > e d , a , Fl,r> w 0 , and obtain from Eq. 
(1?51) the A = 0 rates, 


k d ft k^ L 

b. 

r\j U ~ Ay 


f^+l), 

1 l! r V wo / 

ff(^-l), 

1 Ll V wo / 


(57) 


with negligible left-to-right rate constants. Eqs. (1511) - 
(1551) for the particle current then reduce to 


(Ip°) 

(Ir H ) 


g 2 w0 f _ \ 
r £ r fl Uo 2 )' 
f Am A wg \ 

r L r fl V A/iV ' 


(58) 


Similarly, we derive the particle current noise 

'Ag 4 


(s" u ) = 

(s£ H ) = 

and the Fano factor 


F ho = 


g 2 w o 

T l Tr 

fAg_ 

T l Tr 


Ag 2 


w 


1 - 


- 1 

_w|_ 

A g 4 

A g 2 


F ah = 1 


A/i 2 


w, 


1 . 


2 > 


(59) 


(60) 


Thus, while at low bias and low temperature the DA- 
AH and DA-HO models similarly behave, at high volt¬ 
age fundamental differences are displayed, particularly 
in the current statistics, (see e.g. Fig. [6]). It can be fur¬ 
ther proved that in the DA-AH model higher order cu- 
mulants scale as C n +i/Cn oc 1, while the DA-HO model 
supports Cn+i/C n oc Ag 2 /uj 2 . In contrast, the single im¬ 
purity Anderson Holstein model shows a different scaling 
altogether, C n+ i/C n oc Ag/co^. These three models 
thus display distinct noise characteristics, a useful input 
for identifying the nature of electron-phonon coupling in 
conducting junctions. 
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FIG. 5. (color online) Charge current and differential conductance as a function of voltage bias for the DA-HO model with 
(a)-(b) cuo = 0.1 eV and (c)-(d) coo = 0.05 eV. Parameters are the same as in Fig. [4] but the mode is further coupled to a 
dissipative phonon bath with T p h = 0.05 eV and T p h =Tl = Tr. 


We proceed by presenting simulation results at finite 
temperature and bias, focusing on the large bias limit 
rather than the linear response behavior. We set the 
Fermi energy /t at zero, and adjust the chemical poten¬ 
tials of the leads symmetrically around the Fermi energy, 
Hl = ~HR, A h = Hr ~ hl ■ For simplicity, we align 
the donor and acceptor energies at the same value and 
use £d = £ a = eo — 0.25 eV. The junction is assumed 
symmetric with T = Trr, and the temperature is taken 
rather low, T < uj o,T, eo- We employ g = 0.1 eV for 
the electron-vibration coupling energy; this value may 
seem large given the perturbative nature of our treat¬ 
ment, requiring g/ojQ 1. However, since in the present 
weak-coupling limit the current simply scales as g 2 , our 
simulations below are representative and can be imme¬ 
diately translated to consider other values for g^. Sim¬ 
ulations were performed by evaluating numerically the 
rates, assuming metals with a wide bandwidth D (larger 
than all other energy scales), and an energy-independent 
hybridization T. 

In Fig. [4] we study the DA-AH model and display the 
current and its derivative, the differential conductance, as 
a function of the applied voltage bias. Panels (a) and (b) 
illustrate results with a relatively high mode frequency, 


loo = 0.1 eV, using T = 0.01 eV or 0.1 eV and T=100 
K. We find that when the hybridization energy is small, 
r < w 0 , the current increases in two steps positioned 
at A h ps 2eo and A/i ps 2(eo + wo). These steps are 
clearly resolved as a two-peak structure in the differential 
conductance, see Fig. |U(b). 

The location of these peaks can be reasoned by in¬ 
vestigating the expression for the current, Eq. (1M1) : For 
the given parameters with site energy eo > g and low 
temperatures, the rates k^ L and k^ L dominate the 
current at finite bias whereas k^ R and k^ R are neg¬ 
ligible. As we gradually raise the bias, we find that the 
vibrational relaxation rate k R ^ L significantly increases 
once hr approaches the site energy, fiR = A/i/2 ps eo, as 
the chemical potential precisely sits then within a region 
of a high molecular electronic density of states, reflected 
by the first jump in the current. At further higher bi¬ 
ases, hr ~ e o + wo, the rate k R ^ L is now strongly en¬ 
hanced since an excess energy wo is available for mode 
excitation, producing the second peak. The energy gap 
between the peaks is therefore given by 2wo- At high hy¬ 
bridization energies, F > ojq, the current reaches higher 
values in comparison to the weak hybridization case, and 
it increases monotonically before saturation. However, 
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FIG. 6 . (color online) Fano factor as a function of voltage 

bias for the DA-AH model. Junction’s parameters are same FIG. 7. (color online) Fano factor as a function of applied 
as in Fig. [4j voltage bias for the DA-HO model, with parameters as in 

Fig. m 


in this large T limit the differential conductance reveals 
only a single broad peak centered around A/x ~ 2eo- 

Panels (c)-(d) in Fig. 2] illustrate the behavior of the 
current and the differential conductance when adopting a 
smaller value for the vibrational frequency, coq = 0.05 eV. 
We observe similar features as in the previous larger-wo 
case, only the broadening F now conceals the two-peak 
structure. It is also notable that the first resonance peak 
at A (i = 2eo is higher in magnitude than the second jump 
at A/x = 2(eo + wo). 

Transport in the DA-HO junction is similarly exam¬ 
ined in Fig. [5] The magnitude of the current is higher in 
comparison to the DA-AH case (Fig. [IJ, for both weak 
and strong hybridization energies, due to the availability 
of many additional channels for excitations and relax¬ 
ations, once the bias exceeds the value 2(e 0 +wo). In fact, 
at high bias voltage the charge current diverges and the 
vibrational mode becomes overly-heated, a phenomena 
referred to as “vibrational instability”—. This heating 
effect can be controlled and avoided if dissipation of en¬ 
ergy from the single-molecular vibration to an additional 
bath (besides the metals) is allowe d 37 ’ 41 , see Appendix B 
for details. Mathematically, this divergence is reflected 
by the denominator in the expression for the current: 
At large bias the current-induced excitation rate exceeds 
the relaxation rate (as we break detailed balance at finite 
bias)— To remove the vibrational instability, results in 
Fig. [5] were obtained by attaching the molecular mode to 
a secondary phonon bath with a finite phonon damping 
rate T p h = 0.05 eV and T p h = Tr = Tr. Interestingly, 
inspecting the differential conductance, we observe that 
the peak at A/x = 2 (eo + uio) is more pronounced relative 
to the first peak at A/x = 2eo- This trend is opposite to 
the DA-AH case in Fig. 0] It is explained by noting that 
the second peak corresponds to the opening up of many 
channels for charge transfer in the case of a harmonic 
mode. 

The zero-frequency Fano factor, Eq. m , is inves¬ 
tigated in Figs. [ 6 ] and [7] We find that this mea¬ 
sure strongly reflects the nature of the vibrational mode: 
In the DA-AH model the Fano factor shows a super- 
Poissonian behavior at low-intermediate biases, but in 


the high bias limit A/x > 2(e 0 + wo) it reaches the value 
1, reflecting a Poissonian behavior. In this high bias 
limit, we receive analytically F = —1 + 2, where the first 
(second) term in Eq. (l55l) contributes —1 (2). While we 
cannot offer a fundamental understanding of the involved 
features in Figs. [ 6 j we confirm that they emerge from the 
behavior of the first term in Eq. (15^1) . while the second- 
equilibrium term maintains the value ^2 at the relevant 
range of applied voltage. The DA-HO junction shows a 
very different behavior, see Fig. [7] Here, a Poissonian 
behavior takes place at relatively low biases, A/x < eg, 
but beyond that F is always super-Poissonian, reaching 
high values when many vibrational states participate in 
the conductance. 

Other theoretical studies have confirmed that molec¬ 
ular junctions may reach very high noise levels due to 
electron-vibration scattering processe s 25 ’ 84- — . It should 
be emphasized however that in these calculations large 
F values were materialized under the assumption of a 
strong electron-phonon interaction, while the molecule- 
lead coupling was assumed weak. In contrast, we are 
concerned here with precisely the opposite arrangement: 
weak electron-phonon interaction but arbitrary large 
metal-molecule coupling, and we reach large values for 
F due to the breakdown of the detailed balance relation 
by the applied bias voltage leading to the participation 
of many vibrational states in transpor t 34 ’ 37 . 

To summarize our observations in this Section, the 
current and its noise can reveal information on the vi¬ 
brational mode participating in the transport process, as 
well as provide input on the hybridization strength of 
the molecule to the leads. The differential conductance 
shows a two-peak structure. The separation between the 
peaks corresponds to (twice) the vibrational frequency, 
and the relative peaks’ height can be attributed to mode 
harmonicity. A strong hybridization, T > oxo, smears out 
the double-peak structure to form a single-asymmetric 
feature. Genuine anharmonicity, e.g., in the form of a 
morse potential rather than a two-state system, should 
lead to a differential conductance similar to that obtained 
in the harmonic case, as long as T is greater than the 
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anharmonic energy scale (r > ujq /D e , with D e as the 
dissociation energy in the morse potential). The signif¬ 
icant qualitative differences in the noise characteristics 
between the DA-AH and the DA-HO junctions could as¬ 
sist in identifying the participating “impurity” mode. 


B. Thermoelectric efficiency and its statistics 

1. Large deviation function for efficiency 

In this Section we study the operation of the DA junc¬ 
tion as a thermoelectric engine. We explore the device av¬ 
eraged efficiency under certain conditions and the statis¬ 
tics of efficiency fluctuations, which should play an im¬ 
portant role in small devices as opposed to the bulk. In 
a recent study, Esposito et al. had analyzed the thermo¬ 
electric efficiency statistics in a purely coherent charge 
transport models. Classical models were examined in 
other studies— The DA junction offers a rich oppor¬ 
tunity to examine the thermoelectric efficiency beyond 
linear response, explore the new concept of efficiency fluc¬ 
tuations, and interrogate the role of quantum effects and 
many-body interactions on the operation of a molecular 
thermoelectric engine. The DA junction is particularly 
interesting in this context: As exemplified in Ref.— and 
below in Fig. [HJb) , the macroscopic thermoelectric ef¬ 
ficiency is identical in the DA-AH and DA-HO models; 
only fluctuations of efficiency reveal signatures of molec¬ 
ular anharmonicity. 

To operate the device as a thermoelectric engine we set 
Tl < Tr and fiL > Hr- The macroscopic thermoelectric 
(TE) efficiency qxE is defined as the ratio between the 
averaged power generated by the engine, 

-W = {jjt L -n R ){I p ) (61) 

to the heat absorbed from the hot reservoir, 

Q = (Iq) = (h) - tm(lp)- (62) 

Namely, fj te = . According to the second 

law, the engine’s efficiency is upper bounded, fjxE < q c , 
with q c = 1 — ijr- as the Carnot efficiency. In the language 
of stochastic thermodynamics, corresponding stochastic 
variables can be defined, the results of measurements dur¬ 
ing the time interval t, the fluctuating work — w = —tw 
and input heat flow q = One can further de¬ 

fine the stochastic efficiency for a single realization as 
tlTE = -w/q. 

In our formalism, we obtain the CGF for work and 
heat by going back to the definition of the characteristic 
function, Eq. m■ Rather that using the counting fields 
A p and X e for charge and energy, we make the following 
substitutions, to obtain cumulants for work and heat, 

A e — > Xq, 

A p t XqflR A I^r)- (63) 


A q and X w are conjugate counting parameters for q = 
Hr - HrN r and -w = {n l ~ Hr)N r , respectively. This 
transformation modifies the form of the fluctuation sym¬ 
metry 

Q{X W , Aqr) = Q( — X w + i/?L, —A q + *(/?£ — Pr))i (64) 
which immediately implies (using X w = \ = 0) that— r— 


w 

1 1 

1 

~T~l~ 

VT l 

T R ) q \ 


By invoking the Jensen’s inequality, Eq. (1651) immedi¬ 
ately returns the bound —(w)/(q) < q c , confirming that 
our definitions for q and w are consistent with classi¬ 
cal thermodynamics. In contrast, efficiency fluctuations 
are typically not bounded, and can take arbitrary values 
because of the stochastic nature of small systems. There¬ 
fore, in general, it is useful to investigate the probability 
distribution function to obtain the fluctuating work and 
heat within the interval t, thus he probability distribu¬ 
tion function Pt(q), to observe the value q within t. Ac¬ 
cording to the theory of large deviations, the probability 
function assumes an asymptotic long time form—, 

Pt(v) ~ e~ tJ ^ (66) 

with J{q) identified as the “large deviation function”. 

We rescaled here and below the efficiency by the Carnot 
value, 

V = Vte/Vc ■ (67) 

The upper bound of the efficiency thus corresponds to 
the value q = 1. 

It can be show n 51 ! 53 ' 56 that the large deviation function 
for efficiency can be obtained from Q(X w ,X q ) by setting 
A q = q rj c X it] , and minimizing it with respect to X w , 

J{rf) = -min Q(X W , q q c X w ). (68) 

The CGFs for the AH and HO mode models are given 
in Eqs. (1351) . (1551) . respectively. To study efficiency fluc¬ 
tuations we use the transformation (1631) . and receive the 
LDF from Eq. (l68l) . Note that we do not explicitly eval¬ 
uate the probability distribution function Pt(q). 

It can be proved that J(q) has a single minimum, 
coinciding with the macroscopic efficiency of the en¬ 
gine, and a single maximum, corresponding to the least 
likely efficiency, which equals to the Carnot efficiency, 

2. Gaussian limit: Linear response theory 

In the linear response limit, i.e. close to equilibrium, 
the stochastic work and heat are assumed to be Gaussian 
variables. It is possible then to derive an explicit expres¬ 
sion for the large deviation function, expressed in terms 
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FIG. 8. (color online) (a) Output power P = (ml — g R )(I p ) 
and (b) macroscopic efficiency fj = fjTE/fjc, comparing exact 
results to the value m. Parameters are t a = td = 0.2 eV, 
ojo = 0.01 eV, g = 0.1 eV, T = 0.1 eV, T L = 300 K, T R = 800 
K, r ph = 0. For linear response calculations we define the 
average temperature as T a = ( Tl + T R )/2, the temperature 
difference AT = T R — Tl, and similarly for the chemical po¬ 
tential, g a = (ml + M«)/2 and Am = Mft — Ml- 


of the Onsager’s response coefficients and the thermody¬ 
namic affinities ^ 54 1 56 . The scaled-dimensionless LDF is 
defined as 


J(r,) = j(rj)/S, (69) 

with the entropy production rate S. In the present Gaus¬ 
sian (G) limit it is given b y 52 ’ 54 : 56 


Jc{r]) 


1 (77 + a 2 + ad + adrj) 2 

4 (1 + a 2 + 2ad){rj 2 + a 2 + 2adrj) ’ 


with the dimensionless parameters 


d = 


J yq 


y/LppLqq 



(70) 


(71) 


Here d describes the degree of coupling in the system, a 
is the affinity parameter. Note that in the Gaussian limit 
Jg(v) is bounded, 0 < Jg{v) fs The minimum value 
Jg(t?g) = 0 is obtained at the average (macroscopic) 
efficiency 


VG = 


a{a + d) 


(72) 


(1 + ad) 

The maximum value Jg{Vc = 1) = 1/4 shows up pre¬ 
cisely at the Carnot efficiency, f] c = 1. 

To simulate (170l) we get hold of d and a by extracting 
numerically the coefficients of the linear response-average 
charge ( I p ) and heat currents ( I q ), 


{Ip) — LppAp -\- LpqAq , ( Iq ) — LqpAp T LqqAq. (73) 


Time-reversal symmetry guarantees that L pq = L qp . The 
affinities responsible for the particle and heat fluxes are 
A p = Pl{pr ~ Ml) and A q = /3 L - /3 R , respectively. The 
average entropy production rate, valid in general non¬ 
equilibrium situations, is S = ( I p )A p + ( I q )A c Note 
that we plotted the dimensionless LDF in Figs. l9l and Hill 
when presenting both exact and linear-response results. 
Only in Fig. [Tl]we retract to the unsealed function, when 
comparing different models. 


3. Numerical Results: efficiency statistics 

We investigate numerically the thermoelectric effi¬ 
ciency and its statistics in the DA model, considering the 
effect of mode harmonicity and beyond linear response 
situations. We begin with macroscopic-averaged proper¬ 
ties. Inspecting Eq. (1571) . we note that the numerator 
in the expressions for the average charge and energy cur¬ 
rents are identical in the DA-AH and the DA-HO models. 
Since the macroscopic efficiency is proportional to the ra¬ 
tio of these two currents, we immediately conclude that 
regardless of whether the mode is harmonic/two-state 
system, the same macroscopic efficiency is to be reached, 
at an arbitrary non-equilibrium condition. However, the 
output power takes different values in the two models. 

In Fig. (Sfa) we display the generated power P = 
{hl — Rn){Ip) as a function of bias voltage mm — Rr for 
both DA-AH and DA-HO models. We find that when the 
mode is harmonic the output power can largely exceed 
values reached in a junction with an AH mode due to the 
availability of many additional channels. In Fig. (HKb) we 
examine the efficiency far from equilibrium, and compare 
the exact value to the linear response limit. For the given 
parameters, linear response theory agrees with the exact 
efficiency calculation as long as Am < 0.03 eV. Interest¬ 
ingly, we find that the device can be made more efficient 
in the nonlinear regime, ( fj m 0.75 at ml — M-R ~ 0.14 eV), 
in contrast to the linear response limit (1721) . which is ob¬ 
tained by linearizing the currents around equilibrium, to 
extract the Onsager coefficients. We further recall the 
scaling P oc g 2 , and that fj itself does not depend on g. 

We now turn our attention to the efficiency statistics. 
In Fig. [9] we display the scaled LDF Jah{ii ) for the 
model with an AH mode, (normalized by the entropy 
production rate S). It is obtained from Eq. (17)51) by 
minimizing the analytical form for the CGF (1331) with 
respect to A™. We further compare the exact LDF with 
the Gaussian limit, Jq H (t]), which is obtained from Eq. 
(El by linearizing (numerically) the currents, to obtain 
the parameters d and a , different in general for the DA- 
HO and DA-AH models. J{rj) does not depend on g given 
the normalization with the entropy production rate. 

We examine the efficiency statistics in the three panels 
of Fig. [9]at different-representative values for the applied 
voltage: (a) linear response limit Ag = 0.025 eV, (b) be¬ 
yond linear response, Ag = 0.1 eV, and (c) around the 
maximal value for efficiency, Am = 0.14 eV. We find that 
the minimum value of Jah(v ) corresponds to the macro¬ 
scopic efficiency, and that the Carnot efficiency rj = 1 is 
the least-likely efficiency. We also confirm that Jq H ( g) 
is always bounded between 0 and 1/4, with the upper 
bound reached precisely at the Carnot efficiency 77 = 1. 
In contrast, we find (numerically) that at 77 = 1 the ex¬ 
act LDF satisfies J(l) < 1/4; equality is reached only in 
the linear response limit. When increasing the bias volt¬ 
age the magnitude of efficiency fluctuations grows, and, 
as expected, the Gaussian approximation Jg{v) becomes 
increasingly unreliable. 
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FIG. 9. (color online) Efficiency LDF for the DA-AH model, showing the exact result Jah(v) and the Gaussian limit Jg H {v ) 
at different biases. Parameters are e a = td. = 0.2 eV, u>o = 0.01 eV, g — 0.1 eV, T = 0.05 eV, T p h = 0 ,Tl= 300 K, Tr = 800 K. 
(a) Ap = 0.025 eV, (b) Ap = 0.1 eV, (c) A/x = 0.14 eV. The vertical dashed line identifies the scaled Carnot efficiency r] c =l. 
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FIG. 10. (color online) Efficiency LDF Jho(v) and f° r HO model for different bias voltage. Parameters are same as 

in Fig. © with (a) A/x = 0.025 eV, (b) A/x = 0.1 eV, (c) A/x = 0.14 eV. The vertical dashed line identifies the scaled Carnot 
efficiency x/ c =l. 


In Fig. [10] we display the LDF for the DA-HO model. 
The observed trends are similar to those of Fig. [9] How¬ 
ever, because of the large entropy production rate taking 
place in the DA-HO junction, the normalized J{rf) re¬ 
ceives rather low values. 

We further compare the two models for the vibrational 
mode and plot the unnormalized LDF SJ(rj) in Fig. (fill) . 
As stated before, the macroscopic efficiency coincides in 
these models. However, quite interestingly, the over¬ 
all statistics differs when increasing bias. In the linear 
regime, the DA-HO model performs as an effective two- 
state system at low temperatures and deviations only 
show up at the tail of the distribution. At high bias, the 
DA-AH model suffers more significant efficiency fluctua¬ 
tions relative to the DA-HO model. 

We conclude this section emphasizing central observa¬ 
tions: The statistics of efficiency can reveal information 
on the mode harmonicity, the Gaussian-linear response 
limit becomes highly unreliable at large bias, as expected. 

V. CONCLUSIONS 

We provided a comprehensive analysis of vibrationally 
assisted charge and energy transport in a donor-acceptor 


type molecular junction. Two limiting models were ex¬ 
amined: (a) In the DA-AH junction the vibrational mode 
was highly anharmonic, consisting of a two-level system, 
(b) The mode was taken as harmonic in the DA-HO 
model. Key results are: 

(i) Employing QME and NEGF approaches for the DA- 
AH and the DA-HO models, respectively, we obtained 
analytical expressions for the steady state cumulant gen¬ 
erating functions, Eqs. m and (HkH) . These results are 
valid to second-order in the electron-phonon strength, 
and correct to arbitrary order in the molecule-metal cou¬ 
pling. The CGF furnishes analytical results for charge 
and energy currents in the system, and for fluctuations 
of these quantities. 

(ii) Our analysis establishes that one can reconcile two 
different-eminent quantum transport techniques: QME 
and NEGF. By taking into account scattering processes 
to the same order in perturbation theory, we showed that 
the QME (used here for treating the DA-AH junction) 
and the NEGF approach (DA-HO model) yielded cor¬ 
responding results. Several works had compared trans¬ 
port predictions from these two approaches, showing de¬ 
viations given the different approximations involve d 56 ’ 93 . 
Our work here is unique in demonstrating that one can 
reconcile results from these two techniques by carefully 
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FIG. 11. (color online) The exact LDF SJ(g) (unnormalized) for the DA-AH and the DA-HO models at different biases (a) 
A g = 0.025 eV, (b) A/x = 0.1 eV, (c) A/x = 0.14 eV. Parameters are same as in Fig. l|9|). The vertical dashed line identifies the 
scaled Carnot efficiency g c = 1. 


taking into account corresponding processes. 

(iii) Expressions for the current, output power, as well 
as the Fano factor (noise) were reached, shown to be sen¬ 
sitive to the properties of the vibrational mode. Specifi¬ 
cally, we found that in our model the two-peak structure 
of the differential conductance directly evinces on the 
mode frequency, while the Fano factor definitely reveals 
information on the mode harmonicity/anharmonicity. In 
contrast, the macroscopic (averaged) thermoelectric effi¬ 
ciency was proved to be identical regardless of the mode 
harmonicity, though fluctuations around the averaged 
value were distinct in the two cases. 

(iv) We had employed the DA junction as a ther¬ 
moelectric engine and studied its efficiency fluctuations 
based on the derived CGF for charge and energy trans¬ 
fer. Previous works of efficiency fluctuations were limited 
to classical models^ - — , or to the quantum - purely co¬ 
herent (noninteracting) regime^. Here, in contrast, we 
examined efficiency fluctuations far from equilibrium in a 
quantum many-body model, using a rigorous approach. 

In future studies we will demonstrate the correspon¬ 
dence between the QME an NEGF in the DA-HO model, 
and examine energy harvesting in a double-dot cell with 
three terminals, to examine the quantum photovoltaic 
effect. 
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APPENDIX A: REAL TIME GREEN’S 
FUNCTIONS 


A. Free electron Green’s function 


The free electronic Green’s functions for the leads, with 
the contour times T\ and T 2 , are given as 

9 k{j ll T 2 ) = -i (r c a fe (ri)a^(r 2 )), k G L,R (Al) 

The projection of this contour ordered Green’s function 
to real time generates four different components, namely, 
lesser (<), greater (>), time-ordered (t) and anti-time 
ordered (f) Green’s functions. The lesser and greater 
components are given as (e.g., for l € L), 

-t 2 ) = i{a\{t 2 )ai(ti)) =if L (ei)e Kdt2 ~ tl) , 
gf{h-t 2 ) = -i {ai{ti)a\(t 2 )) = -i[ 1 - /L(ez)]e l£i( ‘ 2_tl) . 

(A2) 


In frequency domain we get 

gf(uj) = 2t Tif L (ei)6(uj - ej), 

9i M = -27t£[1 - f L (ei)\S(u - e/). (A3) 

The following relations between different components of 
the Green’s functions are valid in both time and fre¬ 
quency domain: 

9i = 9i + gf = 9i + gf ) 

g\ = gf-g? = gT-g r n (A4) 

where g\ and g\ are time-ordered and anti time-ordered 
Green’s functions. The retarded Green’s function is de¬ 
fined as 

gl(ti,h) = -i@(h — t 2 )({ai(ti), a] (f 2 )}), (A5) 

and the advanced Green’s function is g1(t\,t 2 ) = 

[gl(t 2 , ti)]*- In a similar manner, counting field depen¬ 
dent Green’s functions are defined on the contour as 


g r {Ti,T 2 ) = -i {T c a r (Ti)al(T 2 )), (A6) 
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where we employ the short notation a r (r) = 
e -i(\ p (T)+e r \ e (T)) j n rea j ^ime, we obtain the lesser 
and greater components, 

gf[ti-t 2 ) = i{al(t 2 )a r (ti)) 

= i fR(e r )e ier< ' t2 ~ tl ' > e l< ' Xp+erXe \ 
g>(ti-t 2 ) = -i(a r (ti)al(t 2 )) 

= -i [1 - /fl(er)]e“ r(ta_tl) e ~ i ( Xp+erXe \ 

(A7) 

In frequency domain they are given by 


< 7 <(w) = 2nifR(e r )8(uj - e r )e l< ' Xp+erK ' > , 

5 >(w) = -2m[l - f R (e r )\S(uj - e r )e~^ Xp+erXe) .(A8) 


Note that relations such as in Eq. (IA41) do not hold for 
counting-field dependent Green’s functions. 


B. Electron-hole Green’s function 

In the main-text we have defined the electron-hole 
propagator in Keldysh space [Eq. (PT3]) ] . In real time we 
receive the four different components as 

F t {t 1 ,t 2 ) = -ig^hfM 2 

l,r 

X [gfih-h) glfo-h) + gKh-h) 

F t (t 1 ,t 2 ) = -ig^hfhrf 

l,r 

X [g\{ti-t 2 )gl{t 2 -ti) + g\{t 2 -ti)gl{tx-t 2 )\, 

F < (t 1 ,t 2 ) = -zg 2 J2\li\ 2 \lr\ 2 

l,r 

X [gf(h-t 2 )g>(t 2 -ti) + g>{t 2 -ti)g<(ti-t 2 )\, 
F>(tl,t 2 ) = -* 5 2 5>| 2 |7r| 2 

l,r 

X [9i{ti-t2)gr{t2-t{] +gf-{t 2 -ti)g>(t 1 -t 2 )\. 

(A9) 

In frequency domain assuming time-translational invari¬ 
ance for the propagator in the steady state limit, these 
components can be written as 


F\u) 

F < (co) 

F>M 



[at (<*>+) Sr (w-) + g\ (w-) gr (w+)], 

dbj r 

~2^ [at ( w +) a ‘ (w-) + g\ (w_) g\. (w+)], 

rh / 

[flf (W+ ) dr (<*>- ) + ^ (W-) g^ (W+ )] , 
[flf (w+) gr (<*>-) + af (w-) g^ (w+)], 


(A10) 


where w± = w' ± Using the relations between the 
Green’s functions [Eq. (1A4[) ] and the expressions for the 


free Green’s functions [Eq. (1A8[) ] we obtain the lesser and 
greater components for the propagator, 


F<(w) = -i2ng 2 ^2 l7l| 2 |7r| 2 /z,(e;)(l - /fl(er)) e ^ x p +t - x ^5{ei - e r - w) 

l,r 

+ ^2 l7i| 2 |7r| 2 /i?(er)(l - fL{ei))e l(Xp+erXe) 8{ei - e r + u) , 

l,r 


(All) 


and F > [oj ) = F < (—ui). Note that at the phonon fre¬ 
quency wo the lesser and greater components of F are re¬ 
lated to the excitation and relaxation rates, defined in the 


QME approach, as E < (w 0 ) = — i k x and F > (oj 0 ) = —i k%- 

The sum of t and t components can be simplified 
following the relations (lA4l) . and using the identity 
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du'g r ’ a (u+)g r ’ a (u}-) = 0, to reach 


/ OO 7 / 

^[.9; < (w+)ff>(w-) +5 ; > ( w -)^(w+) + 5i > (w+)ff<(w_) + ft < (w_)ff>(w + )] 
-oo Z7r 


= F > (w) + F < (w) 


(A12) 


We further obtain 

F*(u;) - F*(w) = 2Re[F r (w)], (A13) 

where 

F r (u j) = ~ig 2 3 4 5 6 ^2,\ll\ 2 \lr\ 2 J duj'^Re[gl(uj + )}g<(u}-) 

l,r 

+ R-e[^(w_)]g ; < (w + )| + 1 -H- r 

(A14) 

APPENDIX B: GENERATING FUNCTION IN 
THE PRESENCE OF A PHONON BATH 

The CGFs derived for the DA-AH and DA-HO models 
can be extended to describe the case with a vibrational 
mode linearly coupled to a dissipative phonon bath. As¬ 
suming this coupling to be weak, it can be shown that the 
formal expressions for the CGFs, Eqs. (l33l) and (1551) . re¬ 
main the same except that the excitation and relaxation 
rates are modified, given now by the sum of electronic- 
baths and phononic-bath induced contributions, 

I.A _ J,A,el i i,pli 
' Si ! 

fcA = feA.el + fc p hj (B1) 

with 

ku = r'ph(^o)'^'ph(^ ? o)5 

k d h = r ph(wo)[l +n p h(u;o)]. (B2) 


^u 'd are the rates defined in the main text, induced 
by the metal leads. Here, r p h(cco) is the coupling en¬ 
ergy of the particular mode u >o to the phonon bath and 
n p h(w) = [exp(/3 p hw) — 1] _1 is the Bose-Einstein distri¬ 
bution function at temperature 1 //3 p h - Note that in the 
presence of the additional phonon bath the fluctuation 
symmetry as written in Eq. (l34l) is not satisfied. To re¬ 
store the symmetry, one should ’count’ as well the energy 
dissipated into the phonon bath. 

For completeness, we include the expressions for the 
charge current and its noise in the dissipative harmonic 
mode model, generalizing Eqs. (1531) and (1551) . 

,jHO\ _ C k uV k d + (kf)'ku , , 

kd — k ’ 

with the short notation (fc® 1 )' = [kfj] R ^ L — [k^] L ^ R and 
(kf)' = [kf]^ 1 " — [kf] L ^ R - The charge current noise is 

(q H 0) = g (^°) 2 , fc® 1 ^ + kfk u + 2 (K'yjkfy 

[ P 1 k d -ku^ k d -ku 

(B4) 


The second term clearly indicates that the noise includes 
terms mixing the effects of the three reservoirs. These 
expressions were used to simulate Figures [5] and [3 taking 
r ph as a constant. 
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